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1.  INTRODUCTION 


This  report  outlines  the  progress  of  our  research  effort  for  1987  and  forecasts  the 
research  plans  for  the  first  half-year  of  activity  of  1988. 

'.Vork  on  the  modeling  of  the  sloshing  fluid  has  resulted  in  a  formulation  and  numerical 
algorithm  for  the  two-dimensional  viscous  sloshing  problem.  On  the  basis  of  our  experiences, 
we  have  chosen  to  use  a  vorticity,  streamfunction  approach  rather  than  a  primitive  variable 
appi  oach  because  of  its  inherent  mass  conservation  property.  A  key  feature  of  the 
formulation  is  the  use  of  a  coordinate  transformation  that  maps  the  fluid  body  into  a  fixed 
geometric  shape.  The  formulation  also  removes  an  initial  singularity  from  the  governing 
equations  that  would  otherwise  cause  the  numerical  method  to  diverge.  Highlights  of  the 
formulation  as  well  as  some  preliminary  results  are  given  in  the  next  chapter. 

We  designed  and  built  an  experimental  test  rig  to  study  the  interaction  between  the 
sloshing  fluid  and  the  spinning  structure.  Results  verified  what  would  be  predicted  from  a 
classical  dynamic  model  of  a  rigid,  torque-free  body.  The  test  rig  exhibited  unstable  motion 
when  spun  about  an  axis  of  intermediate  moment  of  inertia,  but  the  motion  was  stable  for 
rotation  about  an  axis  of  maximum  moment  of  inertia.  Details  are  presented  in  Chapter  3  of 
this  report.  A  dimensional  analysis  was  also  conducted,  to  identify  the  significant  parameters 
associated  with  the  dynamic  behavior  of  such  a  system.  The  results  are  described  in 
Chapter  4. 

On  the  basis  of  experience  with  the  initial  test  rig  configuration,  modifications  were 
incorporated  to  allow  for  rotation  of  the  spinning  assembly  about  an  axis  of  minimum  moment 
of  inertia.  Instrumentation  was  also  added  to  monitor  the  motion  with  a  data  acquisition 
system.  Currently,  the  transducers  are  being  calibrated  and  fine-tuned  for  test  runs 
scheduled  for  the  near  future. 

Development  of  two  versions  of  a  computer  simulation  model  of  the  test  rig  dynamics  is 
also  well  under  way.  Using  equivalent  pendula  to  model  the  sloshing  fluid,  we  will  be  able  to 
make  useful  comparisons  between  predicted  behavior  and  experimental  results.  Then,  when 
the  computational  fluid  model  is  sufficiently  developed,  it  will  replace  the  pendula  in  the 


2.  CFD  MODEL  OF  SLOSHING 


One  of  the  primary  goals  of  this  study  is  to  develop  a  numerical  model  capable  of  describing  the 
sloshing  motion  that  occurs  within  the  fuel  stores  of  a  STAR  48  Communications  Satellite  during  the 
perigee  phase  of  its  orbit.  When  the  power  assist  module  (PAM)  of  the  spacecraft  fires  its  thruster  to 
establish  the  desired  geosynchronous  earth  orbit,  the  fluid  within  the  fuel  stores  experiences  a  sudden 
axial  thu  ust.  Understanding  the  essential  physics  of  the  fluid  reaction  to  this  sudden  thrust  is 
essential  in  establishing  an  accurate  numerical  model.  Consequently,  a  two-dimensional  model  that 
significantly  simplifies  the  geometry  is  used  to  study  the  fluid  motion  that  occurs  in  a  partially  filled 
container  as  a  result  of  the  sudden  application  of  an  external  force.  A  detailed  description  of  this 
model  as  well  as  some  preliminary  results  are  included. 


2.1.  FORMULATION  OFTHE  TWO-DIMENSIONAL  VISCOUS  MODEL 
The  model  consists  of  a  rectangular  container  partially  filled  with  an  incompressible  viscous 
fluid  (see  Fig.  1 ).  Initially,  the  container  is  stationary  and  the  fluid  is  in  a  state  of  hydrostatic 
equilibrium.  At  time  t  =  0  +  ,  an  external  horizontal  force  is  applied  to  the  container  such  that  the 
container  experiences  a  constant  acceleration  toward  the  right.  This  sudden  application  of  the 
external  force  causes  a  singularity  in  the  governing  equations  at  time  t=0  +  .  This  singularity  is  most 
clearly  evident  by  the  instantaneous  adjustment  of  the  pressure  field  (given  the  assumption  of  an 
incompressible  fluid)  to  the  newly  imposed  conditions.  In  order  to  model  correctly  the  fluid  motion 
that  occurs  moments  after  the  application  of  the  external  force,  one  must  remove  this  singularity 
from  the  governing  equations.  This  removal  is  accomplished  by  the  scaling  used  in 
nondimensionalizing  the  governing  equations. 


2.1.1.  The  Governing  Equations 

The  mathematical  rT'ode!  describing  the  fluid  motion,  which  is  induced  as  a  result  of  the  sudden 
container  acceleration,  is  obtained  through  use  of  first  principles  of  fluid  motion.  The  equations 
describing  the  conservation  of  mass  and  momentum  for  the  fluid  are  formulated  for  a  noninertial 
reference  frame  that  moves  along  with  the  container.  These  equations  are  subsequently  used  to 
obtain  the  following  dimensional  equations: 
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Figure  1.  Geometry  of  the  model. 


5 


S 


k 


I 


a  CO  d(f,u)  — 2  — 
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where 


<Xp,q)  dp  dq  dq  dp  —2  d2  a2 

-  -  —  —  _  —  —  and  y  =  -  +  — 

a(x^)  dx  dy  dx  dy  dy 


Equation  (la),  the  vorticity  transport  equation,  describes  *he  conservation  of  vorticity.  Equation  (lb) 
is  the  streamfunction  Poisson  equation,  which  is  obtained  by  expressing  the  curl  of  the  velocity  vector 
in  terms  of  the  streamfunction.  Here,  the  horizontal  and  vertical  components  of  velocity  are  defined 
in  terms  of  the  streamfunction  by 


df  -  df 
a  —  —  —  and  u  =  —= 


respectively.  The  streamfunction-vorticity  approach  was  chosen  over  the  primitive  variable  approach 
because  of  its  inherent  mass  conservation  advantage.  The  pressure  variation  through  the  fluid  is 
obtained  from  the  following  Poisson  equation  that  is  derived  by  use  of  the  continuity  and  momentum 
eq"ations. 


lJx  dy  Xdxby' 

The  motion  of  the  free  surface  is  determined  from  a  kinematic  condition,  which  is  based  upon  the 
assumption  that  fluid  particles  on  the  free  surface  must  remain  attached. 


dR  _  -  -  dR 

dt  dx 


Here,  R(x,t)  denotes  the  free  surface  position  with  respect  to  the  bottom  container  wall. 
Equations  (la,  lb,  and  4)  constitute  the  necessary  set  of  governing  differential  equations  for  the 
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problem.  However,  completion  of  the  model  requires  the  additional  specification  of  initial  and 
boundary  conditions. 


2.1.2.  Initial  and  Boundary  Conditions 


The  initial  conditions  at  t  =  0-  are 


to  =  0  —  f  ,  R(x)  =  constant ,  and  p  —  PQ  —  g  p(R  —  y  1 


(5) 


where  p0  is  the  pressure  at  the  free  surface. 

The  following  boundary  conditions  along  the  container  walls  are  obtained  with  use  of  the  no¬ 
slip  condition  and  momentum  conservation  at  the  walls: 

for  the  streamfunction 


f  =  0  at  x  =0,i  —  a,  and  y  =  0 


(6) 


for  the  vorticity 


and 

for  the  pressure 


du  ?  p 

p  —  =  —  —  —  pg  at  i=0  and  x  —  a 
dx  dy 

d(o  d  p 

p  —  =  —  +  P<?  at  y  =  0 
dy  dx 


dp 

d(j) 

-  Pb  +  p  — 

at  x  =  0  and  x  -  a 

dx 

dy 

if.  - 

3w 

-  P£  -  P  — 

at  y  =  0 

dy 

dx 

(7) 


(8) 


Here,  p  denotes  the  fluid  viscosity  and  p  the  fluid  density.  The  acceleration  due  to  gravity  is  denoted 
by  g,  while  the  acceleration  experienced  by  the  fluid  due  to  the  external  force  application  is  denoted 
by  q.  Note  that  the  vorticity  and  pressure  gradients  are  related  by  the  boundary  conditions  at  the 


7 


solid  walls.  Since  the  flow  is  driven  by  the  sudden  change  in  the  pressure,  it  becomes  essential  that 
the  flow  field  gets  coupled  with  the  pressure  field. 

One  can  obtain  the  boundary  conditions  at  the  free  surface  by  requiring  that  the  normal  and 
tangential  stresses  be  continuous  across  the  liquid  and  vapor  interface.  Since  the  free  surface  is  in 
contact  with  a  much  less  dense  or  viscous  vapor,  one  can  reasonably  assume  that  stresses  on  the  vapor 
side  are  negligible  compared  with  the  ones  on  the  liquid  side.  On  the  basis  of  this  assumption,  the 
following  conditions  have  to  be  satisfied  at  the  free  surface: 


— 

p  —  p  —  2\i  -  (normal  stress) 

0  dn 


(9) 


du  du  _ 

-  +  -  —  k u  =0  (tangential  stress) 

dn  3s  s 


(10) 


where  here  n  denotes  the  normal  and  s  the  tangential  directions  at  the  free  surface.  The  free-surface 
curvature  is  denoted  by  k.  Equation  (9)  is  used  to  determine  the  pressure  variation  along  the  free 
surface,  while  Eq.  (10)  is  used  to  obtain  the  streamfunction  variation  along  the  free  surface. 


2.1.3.  Coordinate  Transformation 

One  of  the  difficulties  associated  with  free-surface  flow  problems  is  that  the  shape  and  position 
of  the  free  surface  are  not  generally  known  a  priori ;  they  are  instead  determined  as  part  of  the 
solution.  This  makes  the  application  of  boundary  conditions  at  the  free  surface  more  difficult  and 
subject  to  inaccuracies  because  of  irregularities  in  the  shape  of  the  surface.  However,  these 
difficulties  can  be  circumvented  if  a  coordinate  transformation  is  introduced  that  immobilizes  the  free 
surface  and  transforms  the  irregularly  shaped  problem  domain  into  a  region  of  a  simple  regular 
shape.  Consequently,  the  following  coordinate  transformation  is  introduced: 

x  =  -  y  =  — r-  .  t  =  t  /V  a/g  ^  ^ 

a  R(x,t) 

which  immobilizes  the  free  surface  and  transforms  the  problem  domain  into  a  unit  square.  The 
coordinate  transformation  causes  the  governing  equations  to  become  more  complex  in  form  as 
information  pertaining  to  the  shape  a^d  motion  of  the  free  surface  is  passed  on  to  the  equations 
through  the  transformation  metrics. 
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2.1.4.  Nondimensionalization 

The  sudden  application  of  the  external  force  acting  on  the  container  at  time,  t  —  0  • ,  is 
transmitted  into  the  fluid  by  means  of  a  pressure  disturbance  that  originates  at  the  container  walls. 
Since  the  fluid  is  being  modeled  as  incompressible,  this  disturbance  is  felt  throughout  the  fluid  at 
time  t  =  0  * .  The  pressure  field  experiences  a  change  in  magnitude  of  order  0(pq)  while  a  weak  fluid 
motion  of  order  0(qt)  occurs  throughout  most  of  the  fluid.  The  core  of  the  fluid  experiences  an 
irrotational,  rigid  body  motion  at  time  t  =  0*;  while  the  fluid  motion  near  the  container  walls  is 
dominated  by  viscous  effects  that  tend  to  resist  the  core  motion  as  fluid  particles  adjacent  to  the  walls 
stick  to  them,  forming  boundary  layers  of  thickness  8.  In  order  to  apply  the  vorticity  boundary 
conditions  accurately,  the  large  gradients  in  the  boundary  layers  must  be  resolved.  To  accomplish 
this,  the  coordinate  normal  to  each  of  the  container  walls  is  stretched  by  a  factor  inversely 
proportional  to  the  thickness  of  the  boundary  layers.  The  following  new  coordinate  transformations 
are  introduced: 

for  the  core 


x  —  2a8 


for  the  left  wall  boundary  layer 


y  —  a8 


x  y  —  a  8 

y=~ 


for  the  right  wall  boundary  layer 


y  -  a  8 


for  the  bottom  wall  boundary  layer 


x  —  2o8  y 

X=~UT’  = 


« .  v 


/  /, 


■mm 


where 


A  =(1-28),  r  =  (B-8),  uP  =  Rt  6=\/n t  IGaV4 ,  Ga  =  — 


g 


Here,  the  superscript  *  indicates  a  stretched  coordinate.  More  generally,  it  will  he  used  to  refer  to  any 
boundary  layer  variable.  The  Galileo  number,  Ga,  is  a  dimensionless  group  that  represents  the  ratio 
of  gravitational  to  viscous  effects. 

In  addition  to  the  three  major  boundary  layer  regions  along  each  of  the  container  walls,  two 
more  regions  are  created  at  the  lower  corners  of  the  container  where  these  boundary  layers  merge. 
These  are  termed  the  left  and  right  corner  layers,  and  the  appropriate  coordinate  waling  for  these 
regions  is  given  by 


for  the  left  corner  layer, 


*  y 

x*  =  —  V*  =  — 

a8  ’  a8 


for  the  right  corner  layer, 


a  —  x  y 

x*  =  -  ,  y*  =  — 

a8  a5 


The  nondimensionalization  of  the  dependent  variables  is  determined  by  scale  analysis.  The 
following  dimensionless  variables  result  for  the  core  and  each  of  the  viscous  layers: 

for  the  core 


3/2  1/2, 
a  g  t 


am«>  P  ~  P0 

«  --  -  p  =  - 

g  mi  P ag 


for  the  boundary  layers 


S'  *7  S'  s'  S'  %  S  S  ■» 


/•  L~’ 


A 


'i  3/2  1/2, 3/2  ”i  1/2, 1/2 

a  g  t  g  t 

for  i  —  l,  r,b,  Ic,  rc 

where  the  subscripts  in  Eq.  (19)  denote  the  following:  l,  a  left  boundary  layer;  r,  a  right  boundary 
layer;  b,  a  bottom  boundary  layer;  Ic,  a  left  corner  layer;  and  rc,  right  corner  layer  quantities.  When 
the  above  coordinate  transformations  and  dimensionless  variables  are  introduced  into  the  governing 
equations,  a  set  of  dimensionless  governing  equations  results  for  each  of  the  six  defined  regions  (see 
Fig.  2).  The  governing  equations  for  the  core  and  the  boundary  layer  along  the  left  wall  are  given 
below: 
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The  governing  equations  for  the  remaining  viscous  regions  are  similar  in  form  to  Eq.  (21).  It 
should  be  noted  that  although  different  scaling  parameters  are  used  to  obtain  the  equations  for 
different  regions,  no  terms  are  neglected  in  any  of  these  equations.  In  addition  to  the  boundary 
conditions  stated  before,  matching  conditions  have  been  developed  for  this  coupled  system  of 
equations.  These  conditions  consist  of  matching  the  value  and  the  normal  gradient  for  each  of  the 
dependent  variables  across  common  boundaries. 


2.1.5.  Dimensionless  Initial  Conditions 


The  initial  conditions  for  each  of  the  six  regions  can  be  obtained  if  the  limit  of  the  governing 
equations  at  t  =  0*  is  considered.  This  yields  the  following  initial  conditions 
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for  the  boundary  layers. 
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where  z*  represents  the  stretched  boundary  layer  coordinate  normal  to  each  wall.  The  equations  for 
the  corner  layers  are  not  given  since  the  solution  to  these  equations  at  t  =  0*  is  trivial.  The  corner 
layers  are  characterized  by  constant  pressure  and  zero  streamfunction  and  vorticity  values  as  no  flow 
occurs  in  the  corner  region  at  t  =  0*  because  of  friction  with  the  walls.  The  above  initial  conditions 
are  subject  to  the  limiting  form  of  the  boundary  and  matching  conditions  as  t-»  0  +  .  The  flow  in  the 
core  is  directly  driven  by  the  pressure  distribution.  This  is  evident  from  the  limiting  form  of  the 
momentum  equations  which  at  time  t  =  0*  become 
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Using  Eq.  (24)  and  the  matching  conditions,  we  can  determine  the  pressure  and  streamfunction 
distribution  in  the  core  analytically  (and  independently  of  the  flow  in  the  boundary  layers).  The  core 
solution  at  t  =  0+  is  given  by 
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where 
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A  similarity  solution  can  be  obtained  for  the  boundary  layers.  The  general  solution  of  the 
boundary  layer  equations  can  be  expressed  in  terms  of  the  complimentary  error  function  as 
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where  C(x,y)  represents  a  constant  for  the  boundary  layer  equations  and  is  determined  by  matching 
with  the  core  along  the  common  boundaries. 

The  initial  conditions  are  also  depicted  graphically  in  Figs.  3,  4,  and  5.  Figure  3  is  a  plot  of 
lines  of  constant  streamfunction.  Note  that  a  clockwise  flow  occurs  from  the  right  to  the  left  container 
wall.  The  closer  the  spacing  between  lines  in  the  vertical  direction,  the  larger  the  magnitude  of  the 
horizontal  velocity  component.  Similarly,  the  closer  the  spacing  between  the  lines  in  the  horizontal 
direction,  the  larger  the  magnitude  of  the  vertical  component  of  the  velocity.  Consequently,  the 
larger  values  of  velocity  occur  near  the  free  surface,  with  the  largest  vertical  component  values 
occurring  near  the  container  walls.  The  value  of  the  horizontal  component  of  velocity  along  the  free 
surface  at  t  =  0  *  is  -qt.  Note,  however,  that  the  horizontal  velocity  component  is  referenced  on  a 
frame  moving  with  the  container;  thus,  in  terms  of  a  stationary  frame  of  reference,  the  stronger  flow 
occurs  near  the  walls.  Figure  4  is  a  plot  of  lines  of  constant  pressure.  Higher  pressure  values  are 
present  in  the  left  half  of  the  container  than  in  the  right  half,  indicating  that  the  free  surface  will  be 
forced  to  rise  along  the  left  half  and  fall  along  the  right  half  of  the  container.  The  similarity  solution 
for  the  flow  in  the  boundary  layers  is  shown  in  Fig.  5.  The  streamfunction,  and  its  first  and  second 
derivatives  with  respect  to  the  boundary  layer  coordinate,  z*,  are  plotted.  In  Fig.  5,  f*z*  is 
proportional  to  the  vertical  velocity  component  in  the  left  boundary  layer,  to  the  negative  of  the 
vertical  velocity  component  in  the  right  boundary  layer,  and  to  the  nega*  ;-t  ~rthe  horizontal  velocity 
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component  in  the  bottom  boundary  layer.  The  vorticity  in  each  of  the  three  boundary  layers  is 
proportional  to  -f*2Z».  Note  that  the  maximum  vorticity  value  occurs  at  the  wall,  and  its  magnitude 
dereases  sharply  to  reach  a  zero  value  at  the  edge  of  the  boundary  layer 


2.2.  The  Sloshing  Model  for  an  Inviscid  Fluid 

In  the  limiting  case,  as  Ga— »<»,  the  fluid  within  the  container  can  be  treated  as  inviscid  This 
greatly  simplifies  the  model  because  in  the  absence  of  any  viscous  effects,  the  need  for  special 
treatment  of  regions  adjacent  to  the  container  walls  no  longer  exists.  The  viscous  boundary 
conditions  at  the  walls  are  replaced  to  allow  for  slip,  and  the  stress  conditions  at  the  free  surface  are 
replaced  by  a  conservation  of  tangential  momentum  requirement.  The  resulting  irrotational  flow  can 
be  described  by  the  solution  of  a  Poisson  equation  for  the  streamfunction,  coupled  with  the  kinematic 
condition  describing  the  free  surface  movement.  The  solution  of  a  Poisson  equation  for  the  pressure 
becomes  optional  because  the  flow  field  is  no  longer  directly  coupled  with  the  pressure  field 

The  fiowfield  obtained  from  this  model  is  not  expected  to  differ  significantly  from  the  one  that 
will  result  from  the  full  viscous  model  for  small  values  of  time.  This  can  be  expected  because  for  small 
values  of  time  the  viscous  effects  are  predominantly  confined  next  to  the  container  walls.  In  the 
absence  of  any  viscous  dissipation,  the  fluid  will  continuously  oscillate  about  the  equilibrium 
position.  However,  as  time  advances,  in  the  presence  of  viscous  dissipation  the  oscillations  will  damp 
out  and  the  fluid  will  come  to  rest  in  a  state  of  hydrostatic  equilibrium  corresponding  to  a  net 
acceleration  due  to  gravity  and  the  external  force.  The  position  of  the  free  surface  at  equilibrium  is 
given  by 


B  (x)  —  B  +Q(J  -x) 

eq  a 

Figures  6  and  7  show  the  transient  behavior  of  the  free  surface  for  the  case  with  B„  =  0.5  and  Q 
=  0.2.  This  corresponds  to  a  square  container  that  is  half  filled  with  a  fluid  experiencing  an  external 
foicc  equivalent  to  0.2  g.  Figure  6  shows  the  position  of  the  free  surface  at  four  different  values  of 
time.  In  Fig.  6a  the  initial  position  of  the  free  surface  is  shown  with  the  solid  line,  while  the  dashed 
line  represents  the  new  equilibrium  position  for  the  free  surface.  At  time  t  =  1.80,  the  inertia  of  the 
fluid  has  carried  it  almost  to  the  maximum  height  that  it  will  reach  along  the  left  container  wall.  The 
fluid  begins  to  descend  along  the  left  wall,  passing  through  the  equilibrium  position  at  approximately 
t  =  2  85,  and  reaches  a  position  close  to  its  initial  distribution  by  t  =  3.90.  Note  that  the  shape  of  the 
free  surface  is  no  longer  fiat.  It  is  suspected  that  this  behavior  is  caused  by  a  nonlinear  interaction 
between  the  various  frequency  components  of  the  disturbance  at  the  walls.  Figure  7  depicts  the 


continuous  departure  of  the  fluid  from  the  equilibrium  position.  A  root  mean  square  measure  of  this 
departure  defined  by 


V  1 2  1 

Drms  =  -  UHx,t)  —  B  (x))~  dx 

Q1  0 

is  plotted  versus  time.  The  peaks  that  appear  in  Fig.  7  represent  the  extreme  departure  positions 
from  equilibrium  that  are  reached  when  the  fluid  rises  to  a  maximum  or  falls  to  a  minimum  height 
along  the  left  wall  The  sharply  shaped  valleys  are  reached  as  the  fluid  passes  through  the 
equilibrium  position  The  sloshing  motion  exhibits  a  period  of  approximately  3.75. 
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3.  TEST  RIG  DESIGN 

One  phase  of  the  project  was  to  design  a  physical  model  to  study  the  effects  of  sloshing  liquid 
fuel  on  the  stability  of  a  spin-stabilized  dynamic  system.  The  design  had  to  allow  for  modification  of 
most  major  dimensions  to  permit  the  testing  of  various  sizes  of  tanks  and  geometric  configurations.  It 
was  decided  that  a  minimum  of  instrumentation  would  be  implemented  until  the  operation  of  the 
system  was  better  understood  so  that  decisions  could  be  made  about  which  parameters  to  instrument. 

A  physical  system  was  conceived  that  consists  of  a  vertical  spin  shaft  with  a  horizontal 
crossbar.  Two  plastic  spheres  are  supported  by  the  horizontal  bar  at  equal  radial  distances  from  the 
vertical  spin  axis.  A  two-degree-of-freedom  Hooke’s  type  universal  joint  is  located  just  below  the 
horizontal  beam  in  the  vertical  shaft  to  allow  the  spin  axis  of  the  horizon.  .  oeam  structure  to  cone.  A 
yoked  sleeve,  hand-actuated  by  a  straight-line-motion,  4-bar  mechanism  is  utilized  to  cover  the 
universal  joint  to  give  initial  stability  and  bending  rigidity  to  the  system  during  spin-up.  A  1/4  hp 
electric  d.c.  motor  is  the  system’s  prime  mover. 

Initially,  we  wanted  to  design  the  system  to  have  a  first-mode  natural  frequency  of  oscillation 
of  0.5  cycles  per  second.  However,  sizing  of  the  spheres  from  the  experimental  results  of  other 
researchers  indicated  that  large  spheres  would  be  required,  including  a  large  mass  of  sloshing  fluid. 

A  compromise  was  achieved  by  an  increase  in  the  natural  frequency  of  the  system  by  a  factor  of  five  to 
2.5  Hz.  This  allowed  the  use  of  readily  available,  6-in. -diameter  spheres  at  a  fill  height  of 
approximately  3.9  in.  The  total  weight  of  one  of  the  sphere  structures,  and  hence  the  beam  end  load, 
is  approximately  five  lb. 

The  sphere-supporting  crossbar  was  analyzed  as  a  non-spinning,  continuous  mass,  cantilever 
beam  with  an  end  load.  Solving  the  frequency  equation  for  the  length  of  beam  resulted  in  a  steel 
beam  27.32  in.  long,  0.625  in.  wide,  and  0.25  in.  thick,  with  a  5.3-lb.  end  mass.  Overall  crossbar 
length  is  54.64  in. 

An  inertial  analysis  was  performed  to  determine  whether  or  not  a  1/4  hp  motor  could  power  the 
system  adequately.  We  decided  to  insert  a  10:1  gear  reducer  after  the  motor  to  utilize  the  higher 
torque  capability  of  the  upper  speed  range  of  the  motor.  Since  the  STAR  48  satellites  are  spin- 
stabilized  at  one  revolution  per  second,  this  was  chosen  as  the  design  speed,  although  speeds  up  to 
three  times  this  are  possible  with  the  motor  and  gear  reducer  system. 

A  right-angle  gear  drive  connects  the  horizontal  output  of  the  gear  reducer  to  the  vertical  shaft 
and  pillow  block  assembly.  The  vertical  shaft  diameter  was  oversized  to  accommodate  pillow  block 
assemblies  that  were  available,  as  well  as  to  allow  for  the  up-scaling  of  future  configurations. 

Figure  8  shows  a  photograph  of  the  completed  assembly  prior  to  any  test  runs.  Note  that  the 
spheres  were  dropped  below  the  center  line  of  the  supporting  beam  in  an  effort  to  locate  the  system 
center  of  mass  at  the  intersection  of  the  two  universal  joint  axes. 


Figure  8.  Photograph  of  completed  test  rig,  initial  design  conf igurat ion 
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Operation  of  the  initial  design  of  the  completed  test  rig  resulted  in  a  very  unstable  motion.  As 
soon  as  the  supporting  sleeve  was  lowered,  the  upper  assembly  would  instantly  hop  over  about  an 
axis  parallel  to  the  supporting  beam. 

An  inertial  analysis  of  the  system  indicated  that  the  spin  axis  was  the  intermediate  axis  of 
inertia.  It  is  well  known  that  spin  of  a  torque-free,  rigid  body  about  its  axis  of  intermediate  moment 
of  inertia  is  an  unstable  motion.  Although  the  test  rig  is  not  rigid,  it  is  assumed  to  be  nearly  so  and 
torque-free  since  the  system  center  of  mass  is  approximately  located  at  the  intersection  of  the 
universal  joint  axes. 

The  differential  equations  for  the  rotational  motion  of  a  rigid  body  can  be  formulated  through 
the  use  of  Lagrange’s  equations.  But  the  first  integrals  of  the  differential  equations  of  motion  can  be 
found  more  directly  in  the  following  manner. 

A  fixed,  inertial  coordinate  system,  XYZ,  located  at  the  system  center  of  mass  was  chosen,  as 
shown  in  Pig.  9,  where  the  support  beam  is  parallel  to  the  X-axis  and  the  Z-axis  is  vertical.  Any  new 
body -fixed  coordinate  system  orientation,  xyz,  relative  to  the  inertial  system,  XYZ,  can  be  reached  by 
a  series  of  three  rotations  about  the  body  axes,  performed  in  the  proper  sequence.  The  Euler  angles  ip, 
0,  and  4>,  which  define  these  three  rotations,  are  presented  here  in  sequence. 

1.  A  positive  rotation  ip  about  the  Z  axis,  resulting  in  the  primed  system  in  Fig.  9. 

2.  A  positive  rotation  0  about  the  y'  axis,  resulting  in  the  double-primed  system  in  Fig.  9. 

3.  A  positive  rotation  <{>  about  the  x*  axis,  resulting  in  the  final  unprimed  system,  xyz,  in 
Fig.  9. 

When  the  orthogonal  projections  of  the  Euler  angle  rate  vectors  onto  each  of  the  axes  of  the  xyz 
coordinate  system  are  determined,  expressions  for  the  body  axis  components  of  the  absolute  rotation 
rate  of  the  body  may  be  written  as 


=  $  —  ip  sin0 
o)^  =  0cos<}>  +  ip  cosQsimj; 
=  ipcos0cos4>  —  0sin4> 


During  initial  startup,  the  test  rig  is  spun  about  a  vertical  axis;  hence,  its  momentum  vector, 
H,  will  be  vertical.  Furthermore,  under  the  assumptions  of  a  torque-free  system,  H  will  remain 
vertical  after  the  supporting  sleeve  is  lowered  and  the  assembly  is  free  to  exhibit  more  complex 
spatial  motion.  If  the  angular  rate  vector  of  ip  is  chosen  to  point  in  a  direction  opposite  the  angular 
momentum,  H ,  the  body-fixed  components  of  the  angular  momentum  arc  found  from  the  projections 
of  H  onto  the  xyz  axes: 


^  **  v.  r  t  L  »  l. 
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All  products  of  inertia  are  zero  for  the  coordinate  system  chosen,  since  two  of  the  three 
coordinate  planes  are  planes  of  symmetry.  Substitution  of  Eqs.  (33)  into  Eqs.  (34)  produces 

—  ipsinG  +  4>  =  (II/I^)sinQ 

tj/cosOsincp  +  Ocos(p  =  (  —  //// ^)cosQsin<p  (35) 

ipcosGcoscJi  —  GsincJ)  =  ( —  II/I  ^)cosQcos*> 

Rearranging  and  solving  for  the  angular  rates,  one  finds 

cu  =  —  II(sin2<b/I  +  cos2<b/I  ) 
yy  zz 

0  =  H(  1//^^  —  l//^)cos0simj)oos4> 

^>  —  FI(V I  —  sin2$/l  —  cos2<$>/I  )sinQ 
xx  yy  22 

These  represent  the  differential  equations  of  motion  in  terms  of  the  Euler  angle  rates. 

Software  was  developed  to  integrate  the  equations  numerically.  A  fourth-order  Runge-Kutta 
routine  was  utilized  as  provided  by  subroutine  RKGS,  a  subroutine  originally  part  of  the  SSP 
subroutine  package  developed  by  IBM.  Angular  rates  were  then  plotted  versus  time. 

To  verify  the  proper  operation  of  the  simulation  software,  we  simulated  the  rigid  body 
dynamics  of  a  STAR  48  satellite  by  utilizing  actual  inertial  characteristics  of  a  STAR  48  after  PAM 
motor  burnout,  where  Ix  =  606,  Iy  =  606,  and  I2  =  318  slug  ft2.  A  constant  spin  of  60  rpm  was 
utilized  and  a  1°  perturbation  was  established  in  the  4>  and  0  directions.  The  results  of  the  simulation 
(Fig.  10)  were  compared  to  actual  telemetered  flight  data  (Fig.  11),  where  the  yaw  gyro  corresponds  to 
<f),  the  pitch  gyro  corresponds  to  0,  and  the  roll  gyro  corresponds  to  ip.  The  frequency  of  the  angular- 
rate  oscillations  for  the  simulation  matched  the  0.5  Hz  frequency  of  oscillation  for  the  STAR  48.  A  90° 
phase  shift  between  the  0  and  4>  angular  rates,  present  in  the  simulation,  indicated  a  coning  motion. 
This  also  agreed  with  the  actual  telemetered  flight  data. 

Simulation  of  the  initial  test  rig  design  was  performed  for  a  five-second  time  interval.  A 
constant  spin  of  60  rpm  was  utilized  and  a  1°  perturbation  was  established  in  the  4>  and  0  directions, 
where  I*  =  0.02,  Iy  =  2.37,  and  Iz  =  2.36  slug  ft2.  Angular  rates  vs  time  are  presented  in  Fig.  12. 
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The  oscillation  o(  the  angu’m  rate  of  4>  about  a  large  negative  value  clearly  indicates  the  instability  of 
this  configuration  and  explains  why  the  test-rig’s  main  supporting  cross  beam  would  instantly  flop 
over  in  the  4>  direction.  This  can  be  visualized  as  the  body  trying  to  reorient  its  spin  from  the  z-axis, 
the  axis  of  intermediate  moment  of  inertia,  to  a  spin  about  the  y-axis,  the  axis  of  maximum  moment  of 
inertia. 

A  computer  routine  was  written  to  analyze  the  inertial  characteristics  of  the  test  rig  for  various 
system  configurations.  Utilizing  this  routine,  we  determined  that  a  shorter  main  support  beam  of 
35.2  in.,  combined  with  the  addition  of  a  54.6  in.  inertial  counter-balance  beam  normal  to  the  main 
beam,  would  result  in  a  maximum  moment  of  inertia  about  the  z-axis  and  hence  a  spin  about  the  axis 
of  maximum  moment  of  inertia.  The  inertias  are  Ix  =  0.13,  Iy  =  2.38,  and  Iz  =  2.46  slug  ft2. 

The  operation  of  the  second  design  configuration  was  simulated  with  the  computer  software 
under  the  same  initial  conditions  as  for  the  first  design.  Angular  rates  versus  time  for  the  second 
design  are  presented  in  Fig.  13.  The  oscillation  of  the  0  and  4>  angular  rates  about  zero  indicates  a 
marginally  stable  motion.  The  simulated  angular  rates  don’t  converge  to  zero  after  they  are 
perturbed  because  the  routine  doesn’t  allow  for  any  energy  dissipation.  The  simulation  results  can  be 
considered  stable  since  the  spin  is  about  the  axis  of  maximum  moment  of  inertia  and  any  real 
physical  system  contains  energy  dissipating  elements  that  would  serve  to  reduce  the  magnitude  of 
the  coning  oscillations  to  zero.  One  must  recall,  though,  that  if  the  spin  were  about  the  axis  of 
minimum  moment  of  inertia,  energy  dissipation  would  serve  to  increase  the  coning  amplitude  until 
the  system  reached  a  fiat  spin-that  is,  spin  about  the  axis  of  maximum  moment  of  inertia. 

Operation  of  the  second  configuration,  shown  in  Fig.  14,  proved  to  be  very  stable.  Any 
perturbation  disappeared  within  one  or  two  seconds.  It  was  impossible  to  determine  whether  any 
fluid  slosh  was  present  and  whether  there  was  any  interaction  between  the  fluid  and  the  physical 
structure. 

It  was  decided  that  a  spin  about  the  minimum  axis  of  inertia,  a  spin  that  is  stable  only  if  there 
is  no  internal  energy  dissipation,  would  best  demonstrate  the  destabilizing  effects  of  sloshing  liquid 
fuel  on  a  spin-stabilized  vehicle. 

Because  of  the  physical  constraints  of  the  test  rig,  no  configuration  with  the  spin  axis  as  the 
axis  of  minimum  moment  of  inertia,  and  with  the  system  center  of  mass  located  at  the  intersection  of 
the  universal  joint  axes,  was  possible.  It  was  determined  that  redesign  of  the  test  rig  would  be 
required  to  thoroughly  study  a  torque-free  spin  about  the  minimum  axis  of  inertia. 


psidcjt,  deg/sec  I  thetaoqt,  deg /sec  I  phidot,  deg /sec 


Figure  14.  Photograph  of  completed  te.st  rig,  second  design  configuration 
spin  about  axis  of  maximum  moment  of  inertia. 
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4  DIMENSIONAL  ANALYSIS 
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In  order  to  maximize  the  usefulness  of  the  test  seL-up,  principles  of  dimensional  analysis  and 
similitude  were  applied  to  the  test  rig.  The  results  of  this  effort  can  be  used  to  design  future  test  rig 
configurations  and  to  scale  test  results  to  obtain  predictions  of  actual  STAR  48  satellite  motions. 

A  list  of  parameters  was  compiled  to  describe  the  system  dynamics.  The  parameters  and  their 
respective  fundamental  dimensions  follow: 
h  =  fluid  depth,  L 

a  =  spin  to  transverse  moment  of  inertia,  assuming  axisymmetric  bodies  only,  l^/It 

t  =  time,T 

on  =  initial  rate  of  spin,  T1 

0  =  nutation  angle 

p  =  fluid  density,  ML*3 

p  =  fluid  viscosity,  ML1T1 

R  =  tank  radius,  L 

g  =  gravitational  and/or  inertial  acceleration,  LT  - 
s  =  liquid  surface  tension,  MT-2 

u  =  displacement  of  sphere  due  to  transverse  beam  flex,  L 
k  =  transverse  spring  rate  of  support  beam,  MT  2 
c  =  effective  beam  damping  in  transverse  direction,  MT  1 
M  -  ratio  of  liquid  mass  to  rotating  mass 
F  =  slosh  force,  response  variable,  MLT  2 
where  M  =  mass,  L  =  length,  and  T  —  time 

Inspecting  the  list  of  parameters,  one  finds  that  many  of  the  terms  are  already  dimensionless. 
These  terms  may  be  considered  u-terms  to  be  used  for  scaling  purposes.  The  first  of  these  terms  are 


nl=0 


n2=0 
n3  =  Af 


All  of  the  remaining  parameters  are  described  by  the  three  fundamental  dimensions-mass, 
length,  and  time  The-  .fore,  the  list  of  parameters  can  be  reduced  by  three  through  the  use  of  the  "pi 
theorem  ” 
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[because  n  terms  are  products  or  quotients  of  the  original  parameters,  and  these  products  or 
quotients  are  dimensionless,  a  general  statement  of  dimensional  homogeneity  of  the  original 
parameters  can  be  formed: 


r^jl  a2  a3  a4 rJa5  ,a6  a7 Lati  a9LaW  ul  1  ,al2  d  , , 0  ,  0  ,,,0 
rcppRgskuh  to  t  =M  L  I 


(38) 


where  the  exponential  a’s  are  integer  constants  and  the  equals  sign  with  a  "d”  over  it  means 
"dimensionally  equivalent.”  Step-by-step  application  of  the  method  yields  the  following 
dimensionless  n-terms: 


II  =  (0/ 

4 


n5  =  F/pftV 


ng  =  c/p  R  to 


n?  =  p/p/?2w 


rig  =  g/R  to 


ng  =  s/p/?3to2 


nw=k/pR3J 


iiu  =  uJR 


(39) 


rl12  =  hJR 


Many  of  these  terms  are  identical  to  those  obtained  by  other  researchers  of  sloshing  fluids.  By 
utilizing  properties  of  n-term  manipulation,  we  can  modify  some  of  the  other  terms  and  find  them  to 
be  equivalent  to  the  results  of  others.  For  example,  multiplying  117  by  (ii#)  1  yields 


,  ,-l  _  H 

n7  M8  0  3/2  1/2 


(40) 


This  term  is  commonly  referred  to  as  the  Galileo  number.  It  represents  the  ratio  of  viscous  effects  to 
gravitational  effects.  For  proper  force  scaling  the  Galileo  number  should  be  correctly  scaled. 

114  is  known  as  dimensionless  time. 

115  is  the  dimensionless  term  relating  the  force  response  parameter,  F,  to  the  inertial  force  of 
the  liquid  mass  This  term  can  be  interpreted  as  a  magnification  factor  for  the  fluid  slosh  force 
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The  beam  spring  rate  and  damping  ratio  were  included  in  the  dimensional  analysis  for  future 
use  when  the  interaction  of  the  elastic  structure  and  sloshing  liquid  is  considered,  n-terms  1 0  and  6 
should  be  used  to  scale  beam  spring  rate  and  damping  ratio. 

117  is  commonly  known  as  Stoke’s  number,  relating  the  viscous  effects  to  inertial  effects. 

(rig)-1  is  known  as  the  Froude  number  and  is  the  primary  scaling  factor  for  the  test  rig.  It 
indicates  the  ratio  of  inertial  to  gravitational  effects.  Hence,  it  indicates  the  way  that  gravity  should 
be  scaled.  Since  all  frequencies  must  scale  as  dictated  by  the  Froude  number,  we  can  scale  the 
prototype  natural  frequency,  (fn)p,  from  the  measured  model  natural  frequency,  (fn)m,  by  requiring 
the  Froude  number  to  be  the  same  for  the  model  and  prototype: 


R(fr 


R(fJ 


(41) 


Rearranging  and  solving  for  the  prototype  natural  frequency,  one  finds 


(f  ) 

'  n  p 


(42) 


Hence,  a  direct  scaling  relationship  has  been  determined  that  relates  the  frequencies,  tank  radii,  and 
effective  gravitational  accelerations  between  the  model  and  prototype. 

Equation  (42)  may  be  used  to  scale  frequencies  for  two  modes  of  free  surface  fluid  slosh.  One 
mode  corresponds  to  oscillations  in  the  vertical  or  ti  ansverse  plane,  while  the  other  mode  corresponds 
to  oscillations  in  the  horizontal  or  circumferential  plane. 

For  the  STAR  48  prototype),  the  effective  gravitational  acceleration  is  simply  the  normal 
component  of  angular  acceleration  since  the  acceleration  due  to  gravity  is  negligible  in  orbit.  For  the 
test  rig  (model),  the  effective  gravitational  acceleration  is  approximated  by  the  vector  combination  cf 
the  normal  component  of  angular  acceleration  and  the  acceleration  due  to  gravity.  Since  the  two 
components  of  acceleration  are  perpendicular  to  each  other,  I  he  magnitude  of  the  effective 
gravitational  acceleration  for  the  test  rig  can  be  determined  in  the  following  form: 


*1/2 


*m  = 


r  (w  ) 


+  (g  r 


(43) 


where  rm  is  a  characteristic  length  that  depends  upon  the  mode  of  oscillation  as  well  as  the  position  of 
the  tank  center  and  its  liquid  contents,  and  g0  is  the  local  gravitational  acceleration  acting  in  the 
plane  of  oscillation.  If  we  utilize  Eqs.  (42)  and  (43)  it  is  possible  to  scale  test-rig  spin  rate,  tank  radii, 
and  slosh  frequencies  to  or  from  the  STAR  4d  spin  rate,  tank  radii,  and  slosh  frequencies: 


V  = 
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[r(w  )n/f) 

P  P  m 


(R)  r  (co  )2  +  (g  f 

p\  mm  ^  o 


*  V 


where  rp  is  defined  for  the  prototype  in  the  same  manner  as  rm  is  defined  for  the  model. 

Consideration  of  the  other  n-terms  reveals  that  multiplication  of  iiq  by  (ny)-'  yields 


This  term  is  known  as  the  Bond  number.  It  indicates  the  ratio  of  surface  tension  forces  to 
gravitational  forces.  Generally,  the  gravitational  forces  are  so  large  compared  to  the  surface  tension 
forces  that  the  Bond  number  need  not  be  scaled.  Only  when  g  is  very  small,  as  it  is  for  a  low-gravity 
space  maneuver  of  a  non-spinning  vehicle,  is  the  Bond  number  an  important  parameter. 

n-terms  1 1  and  12  are  geometric  scaling  relationships.  Term  12  represents  the  ratio  of  tank  fill 
height  to  tank  radius  and  is  an  integral  part  of  the  calculation  of  fluid  slosh  resonances. 

Theoretically,  total  similarity  between  the  model  and  prototype  can  only  be  accomplished  by 
satisfying  all  of  the  dimensionless  n-terms  between  the  model  and  prototype.  Unfortunately,  this  can 
be  impossible  for  some  of  the  terms  because  of  various  physical  constraints.  For  example,  the  Galileo 
number  is  very  difficult  to  satisfy  because  of  the  small  viscosities  of  most  propellants  and  the  large 
tank  radii  of  most  prototypes.  Therefore,  the  modeling  of  a  large  tank  with  a  smaller  tank  could 
require  the  use  of  an  exotic  or  nonexistent  fluid.  Consequently,  researchers  must  decide  which  of  the 
terms  to  satisfy  in  order  to  develop  an  accurate  model  of  the  desired  physical  parameters. 


5.  CURRENT  STATUS  AND  FUTURE  GOALS 


Presently,  work  is  continuing  in  three  specific  directions.  Development  of  the  fluid  model 
remains  a  high-priority  effort;  redesign  and  instrumentation  of  the  test  rig  is  nearing  completion;  and 
computer  simulation  of  the  rotating  structure,  with  the  fluid  modeled  as  equivalent  pendula,  is 
progressing. 

The  work  on  the  two-dimensional  rectangular  model  is  expected  to  be  completed  by  the  end  of 
January  1988.  By  that  time  more  data  will  have  been  gathered  from  the  inviscid  model.  Debugging 
work  on  the  viscous  model  should  be  complete  and  the  model  expected  to  yield  useful  results  in  the 
near  future.  Work  on  the  three-dimensional  general  model  will  intensify  in  early  February,  and  we 
expect  that  useful  results  can  be  obtained  by  early  summer. 

The  test  rig  has  been  redesigned  to  raise  the  universal  joint  higher  above  the  supporting  frame. 
This  will  allow  the  newly  modified  structure  to  be  rotated  about  its  axis  of  minimum  moment  of 
inertia.  Transducers  have  been  obtained  to  monitor  the  spin  rate  of  the  main  vertical  shaft,  the  two- 
dimensional  rotation  of  the  structure  with  respect  to  the  universal  joint,  the  sloshing  of  the  liquid  in 
the  spherical  tanks,  and  the  elastic  vibration  of  the  rotating  structure.  All  transducers  have  been 
mounted  on  the  test  rig  except  for  the  strain  gages,  which  are  to  be  used  to  observe  the  vibrational 
motion. 

The  input  speed  is  being  measured  with  a  DC  generator  mounted  on  the  speed  reducer  between 
the  motor  and  main  vertical  shaft.  Two  potentiometers  are  used  to  follow  the  spatial  motion  of  the 
rotating  assembly  relative  to  the  universal  joint,  and  three  photo-potentiometers  are  positioned  on 
each  of  the  transparent  spheres  to  track  the  location  of  the  liquid  surface  during  test  runs.  A  set  of 
slip  rings  transfers  signals  from  the  rotating  member  to  the  computer  acquisition  system  used  to 
record  data. 

Calibration  curves  for  the  various  transducers  are  presently  being  developed,  and  different 
lighting  schemes  to  improve  the  output  from  the  photo-potentiometers  are  being  evaluated.  The 
strain  gages  will  be  added  after  the  initial  set  of  runs  with  existing  instrumentation.  The  system 
should  be  debugged  and  reauy  to  produce  usable  experimental  results  by  the  beginning  of  February. 
Its  present  status  is  depicted  in  Fig.  15. 

Meanwhile,  two  computer  models  are  being  developed  for  simulation  of  the  test  rig.  A  dynamic 
analysis  software  package  is  available  locally  that  will  allow  the  test  rig  assembly  with  the  fluid 
replaced  by  two-dimensional  pendula  to  be  described  for  dynamic  analysis  without  our  having  to 
derive  the  actual  equations  of  motion.  Results  from  this  model  should  be  available  by  mid-February. 
In  addition,  the  equations  of  motion  are  being  developed  by  use  of  a  Lagrangian  approach.  These  will 
then  be  solved  by  one  of  the  differential  equation  subroutines  readily  available  today.  Results  from 
the  two  computer  models  will  be  compared  with  each  other  for  verification  purposes,  and  they  will 
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also  be  compared  with  experimental  results.  The  equivalent  pendula  will  eventually  be  replaced  by 
the  three-dimensional  fluid  model  for  simulation  purposes. 
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